The Quasi-Reversibility Method for the 
Thermoacoustic Tomography and a 
Coefficient Inverse Problem 

Michael V. Klibanov*, Andrey V. Kuzhuget*, 
Sergey I. Kabanikhin**, and Dmitriy V. Nechaev v 
^Department of Mathematics and Statistics 
University of North Carolina at Charlotte, 
Charlotte, NC 28223, USA 
** Sobolev Institute of Mathematics 

of the Siberian Branch 
of the Russian Academy of Science 
Prospect Acad. Koptyuga 2, 
Novosibirsk, 630090, Russia 
v Lavrent'ev Institute of Hydrodynamics 
of the Siberian Branch 
of the Russian Academy of Science 
Prospect Acad. Lavrent'eva 15 
Novosibirsk, 63090, Russia 
E-mails: mklibanv@uncc.edu; akuzhuge@uncc.edu; 
kabanikh@math.nsc.ru; nechaev@hydro.nsc.ru 

February 2, 2008 



Abstract 

An inverse problem of the determination of an initial condition in a hyperbolic 
equation from the lateral Cauchy data is considered. This problem has applications 
to the thermoacoustic tomography, as well as to linearized coefficient inverse problems 
of acoustics and electromagnetics. A new version of the quasi-reversibility method is 
described. This version requires a new Lipschitz stability estimate, which is obtained 
via the Carleman estimate. Numerical results are presented. 
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1 Introduction 

In this paper we propose a new version of the Quasi-Reversibility Method (QRM) for the 
inverse problem of the determination of an initial condition in a hyperbolic equation from 
the lateral Cauchy data. We discuss applications of this inverse problem to thermoacoustic 
tomography, as well as to linearized coefficient inverse problems of acoustics and electromag- 
netics. Using the Carleman estimate, we prove convergence of our version of the QRM. We 
also present numerical results. In particular, we show that this version of the QRM enables 
one to image 5— like functions, i.e., narrow high peaks. 

Let f2 C W 1 be a convex domain with a piecewise smooth boundary <9f2 and 2R be the 
diameter of Q,2R = max. xye fi \x — y\ . Let T = const. > R. Denote Qt = fix (0, T) . 
Consider the elliptic operator L(x,t) of the form 

n 

L(x, t)u = Au + bj (x, t) Uj + b (x, t)u t + c (x, t) u, 

3=1 

where Uj := d Xj u. We assume that all coefficients of the operator L belong to C (Q T ) ■ Let 
the function u G H 2 (Qt) be a solution of the hyperbolic equation in the cylinder Q T , 

u tt — L(x,t)u + F(x,t) in Q T , (1.1) 

F G L 2 (Qt) with initial conditions 

u (x, 0) = (x) , ut (x, 0) = i\) (x) , <p e H 1 (Q) , i\) e L 2 (Q) . (1.2) 

We consider the following 

Inverse Problem 1. Let one of functions (p or ip be known and another one be unknown. 
Determine that unknown function assuming that the following functions / and g are given 

du 

u\ ST =f(x,t), — \s T =g(x,t), S T = dnx(0,T), (1.3) 

where v is the unit outward normal vector at dQ. We call the problem of the determination 
of the function tp the "^—problem" and the problem of the determination of the function ip 
the "■?/>— problem" . 

In principle, in the case T > 2R one should not assume that one of functions (f or ip is 
known. This is because for T > 2R the following Lipschitz stability estimate takes place (see 
[4], [10], [11] and Theorem 2.4.1 in [13]) 

II«IIhi ( o t ) < c (||/|| h i (St) + \\g\\ L2{ s T) + II^IL^q^) • 
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Here and below C denotes different positive constants depending only on Q, T and C (Q T ) 
norms of coefficients of the operator L. However, since numerical studies for the case of 
the finite domain were conducted in previous publications [4], [12], we are interested here in 
solving the Inverse Problem 1 in an unbounded domain, which was not done before. This 
leads us to the case T > R. Namely, we want to solve an analogue of the Inverse Problem 1 in 
a quadrant, assuming that the lateral Cauchy data are given only on parts of two coordinate 
axis. We are motivated by the publication [14] , where the Lipschitz stability was proven for 
an analogue of Inverse Problem 1 for the case of either a quadrant in M 2 or a an octant in IR 3 , 
assuming that one of initial conditions (1.2) is zero, and the second one has a finite support. 

We now specify conditions of our numerical study. Suppose that equation (1.1) is ho- 
mogeneous with F (x, t) = and it is satisfied in D T = IR 2 x (0, T) . Consider the quadrant 
QU = {xi,X2 > 0} . And also consider the square SQ C QU, 

SQ (a) = {0 < xi, X2 < a} . 

Suppose that 

<p(x) — if) (x) = outside of SQ (a) . (1.5) 
Then the energy estimate implies that 

u (x, t) = 0, V (x, t) G {x | x G QU, dist (x, SQ (a)) > T} x (0, T) . (1.6) 

Denote 

Tit = fa G (0, a + T) , x 2 = 0} x (0, T) , 

r 2 r = {x 2 e (0, a + T) , x x = 0} x (0, T) , 
T 3 t = {xi = a + T, x 2 G (0, a + T)} x (0, T) , 
T 4T = {x 2 = a + T, X! e (0, a + T)} x (0, T) , 
see Figure [H Then by (1.6) 

du 

u = — — = on T 3 t U Tzit- (1-7) 
av 

Hence, we focus our numerical study on 

Inverse Problem 2. Let equation (1.1) be satisfied in D T with initial conditions (1.2) 
satisfying (1.4). In this case Q := SQ (a + T) .Suppose that one of these initial conditions 
is zero. Determine the second initial condition, assuming that functions / and g are known, 
where 

du 

u |r 1T ur 2T = / (x,t) , — |r 1T ur 2T = g {x,t) . (1.8) 

Suppose for a moment that only the function / (x, t) is given. Then one can solve the 
boundary value problem for equation (1.1) with F = outside of the square SQ(a + T) 
with the following initial and boundary data 

u (x, 0) = u t (x, 0) = 0, x G R 2 \SQ(a + T), 
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Figure 1: Geometry for the Inverse problem 2. 



u |ri T ur 2 T— / ( x i t) , u |r 3T ur 4 T— 0- 
This gives one in a stable way the normal derivative g (x, t) on Fit U T 2 t- Thus, we arrive 



at Inverse Problem 2. It was proven in [14] that if 



T > 



i\/2 



V2 



1.9) 



and one of functions ip or tp equals zero, then the following Lipschitz stability estimate is 
valid 



\u\ 



HHG T ) 



HHT T ) 



+ \\g\ 



L 2 (r T ) 



)■ 



(1.10) 



where Tn 



T = r 1T UT 2T and ||/|| H i (rT) = II/IIhi^) + 11/11^^) ■ The estimate (1.10) implies 
a similar estimate for the unknown initial condition [14]. The knowledge of the fact that 
one of initial conditions was zero was used in [14] for either odd or even extension with 
respect to t of the function u(x,t) in {t < 0} , depending on which of initial conditions 
was assumed to be unknown. The proof of [14] is based on the Carleman estimate. The 
method of Carleman estimates was first applied in [11] to obtain the Lipschitz stability for 
the hyperbolic problem with the lateral Cauchy data, also see [10] and Theorem 2.4.1 in [13]. 
Prior to [11] the Lipschitz stability for the hyperbolic problem with the lateral Cauchy data 
was obtained in [20] by the method of multipliers, but only for the case when lower order 
terms in (1.1) are absent. The use of the Carleman estimate enables one to incorporate lower 
order terms and also to extend to the case of hyperbolic inequalities. Recently the method 
of [10], [11], [13] was extended to hyperbolic equations with the non-constant principal part, 
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see, e.g. [22]- [24]. The method of multipliers was recently extended to the case of non-zero 
lower order terms in Theorem 3.5 of the book [6]. 

The above problems were previously solved numerically in [4], [7], [12] and [15]. The work 
[7] was the first one, where the problem of thermoacoustic tomography was formulated and 
solved numerically as Inverse Problem 1, i.e., as the hyperbolic Cauchy problem with the 
lateral data. The QRM for the latter problem was used in [4]. The QRM was first proposed 
in the book [19] for a variety of ill-posed boundary value problems. Its convergence rates 
were established in [9] and [11] for the cases of Laplace and hyperbolic equations respectively 
and in section 2.5 of [13] for elliptic, parabolic and hyperbolic equations. In particular, it 
was shown in [13] that one can work with weak H 2 solutions of QRM instead of strong H 4 
solutions of the original book [19]. Also, see [2] and [3] for the recent results for the QRM 
for the elliptic case and [8] for the application of the QRM to linearized coefficient inverse 
problems for parabolic equations. The main tool of works [4], [8], [9], [11] and [13] is the 
tool of Carleman estimates. 

There are three main differences between the current paper and the previous works on the 
QRM for hyperbolic equations. First, we take into account boundary conditions via including 
them in the Tikhonov regularizing functional J e . Unlike this, boundary conditions were made 
zero in [4] via subtracting off a corresponding function, and they were treated via integration 
by parts in [12]. Second, we incorporate in J £ a penalizing term, which reflects our knowledge 
of one of initial conditions. We show numerically that without this term we cannot image 
well maximal values of the unknown initial condition inside of narrow peaks. On the other 
hand, since cancerous tumors can be modeled as narrow peaks, it is interesting to image 
those peaks in the application to thermoacoustic tomography considered below. These first 
two ideas for J e are taken from [15]. Mainly because of the second difference we cannot 
apply previously derived convergence results and thus, need to prove convergence of our new 
version of the QRM. In particular, we need to prove a new Lipschitz stability estimate 
(Theorem 4.1). Finally, the third difference is that while H 2 finite elements were used in 
[4] and [12], we use finite differences now. Note that while smooth slowly varying functions 
were reconstructed numerically in [4] and [12], our numerical experiments reconstruct both 
those functions and «5 — like functions. «5 — like functions were also reconstructed in [15] for 
the Inverse Problem 2. However, the numerical technique of [15] is different from one of the 
current paper. The method of [7] and [15] is based on the representation of the function 
u(x,t) via truncated Fourier series and minimization of the resulting residual least squares 
functional. 

In section 2 we describe applications of Inverse Problems 1 and 2. In section 3 we describe 
the version of the QRM we use here. In section 4 we prove a new Lipschitz stability estimate. 
In section 5 we prove convergence of our method, based on the result of section 4. In section 
6 we describe our numerical implementation. In section 7 numerical results are presented. 
Conclusions are drawn in section 8. 
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2 Applications 



In this section we discuss two applications of above inverse problems 

2.1 Thermoacoustic tomography 

Inverse Problems 1 and 2 arise in thermoacoustic tomography [4], [16], [25]. In this case 
the target is subjected to a short electromagnetic impulse. The electromagnetic energy is 
absorbed. As a result, temperature is increased and the target is expanded. This causes 
a pressure wave, which is measured as a change in the acoustic field at the boundary of 
the sample. Assuming that the absorption of the electromagnetic energy is spatially varying 
inside the sample, the resulting wave field is carrying the signature of the inhomogeneity. On 
the other hand, cancerous regions absorb more than surroundings. This leads to applications 
in medical imaging. Hence, the problem is to calculate the absorption coefficient a(x) of 
the sample using time dependent measurements at its boundary. Let (3 be the thermal 
expansion coefficient, c p be the specific capacity of the medium and Iq be the power of the 
source. Usually f3, c p and Iq are known. Also, assume that the speed of sound in the medium 
is constant and equals 1. Let u(x, t) be the pressure wave. It was shown in e.g., [4] that 

u tt = Au,(x,t) GR 3 x (0,T), (2.1) 
u(x,0) = a(x)I —, u t (x,0)=0. (2.2) 

Cp 

Suppose that we measure the function u(x, t) at the boundary of the domain Q and a(x) =0 
outside of fl. Then those measurements give us the boundary value problem for equation 

(2.1) outside of Q with zero initial conditions. Solving this problem, we uniquely determine 
the normal derivative of the function u(x,t) at dfl. Thus, we arrive at the tp— problem. 

A different approach to the problem of thermoacoustic tomography is currently actively 
developed in a number of publications. In this approach the solution of the problem (2.1), 

(2.2) is presented via the Poisson-Kirchhoff formula, which leads to the problem of integral 
geometry of recovering a function via its integrals over certain spheres, whose centers run over 
a surface and radii vary. Then uniqueness theorems are proven for this case and inversion 
formulas are derived, see, e,g., [1], [5], [16], and [17]. In particular, works [1] and [16] include 
the case of a variable speed and [16] and [17] include numerical examples. A survey of these 
developments can be found in [16]. Also, see §1 of Chapter 6 of the book [18] for an example 
of the ill-posedness of this integral geometry problem for the case when centers of spheres 
run over a plane in R 3 . 

2.2 Linearized inverse acoustic and electromagnetic problems 

There is also another application, in which Inverse Problems 1 and 2 can be considered as 
linearized inverse acoustic and inverse electromagnetic problems. The idea of this subsection 
is motivated by §1 of Chapter 7 of [18] and §3 of Chapter 2 of [21]. We present this application 
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now without discussing delicate details about the validity of the linearization. In this setting 
the point source is running along a surface and time dependent measurements of back- 
reflected data are performed at the positions of the source. In [11] the Newton-Kantorovich 
method was presented for the case when the source position is fixed and the time dependent 
measurements are performed at a surface. 

Let the function ct(x) G C (R 3 ) be strictly positive, a (x) > const. > 0. Consider the 



Cauchy problem for the hyperbolic equation 

a (x) W U = A X W + 47T<5 {x — Xq, t),(x,t) GR 3 x (0,T), (2.3) 

w(x,x o ,0) = w t (x,x ,0) = 0, (2.4) 

where xq G R 3 is the source position. It is well known that in acoustics a (x) = c~ 2 (x) , 



where c (x) is the speed of sound in the medium, and in some situations of the electromag- 
netics a (x) = (fie) (x) , where /i and e are respectively magnetic permeability and electric 
permittivity of the medium. We pose the following 

Inverse Problem 3. Suppose that the function a (x) — 1 outside of the domain f2 and 
it is unknown inside of this domain. Determine this function for x G f2, assuming that the 
following function p (xo, t) is known 

w (x , x , t) \ xoe aa= V (^o, *) • (2.4) 

The full Inverse Problem 3 is difficult to address because of its nonlinearity. Hence, we 
consider now a linearized problem. Similarly with §3 of Chapter 2 of [21], suppose that the 
function a (x) can be represented in the form 

a (x) = 1 — £a (x) , 

where £ G (0, 1) is a small parameter. Hence, the term £a (x) is a small perturbation of 1. 
We assume that this perturbation is unknown, i.e., the function a (x) is unknown. Again, 
similarly with [21], we can formally set at £ — > 

w (x, x , t) = w (x, x , t) + £u>i (x, x , t) + (£ 2 ) , (2.5) 

where functions wo and w\ are independent on £. This setting was rigorously justified in §3 
of Chapter 2 of [21] for the case of the telegraph equation 

w tt = Aw + (a (x) + £ai (x)) w. (2.6) 

It was also justified in §1 of Chapter 7 of [18] for equation (2.6) without the introduction of 
the parameter £, which is actually introduced here for convenience only. Indeed, instead, 
one can assume that a (x) — 1 — a (x) , where \a (x)\ « 1. 

Substituting (2.6) in (2.3) and (2.4) and dropping the term with O (£ 2 ) , we obtain that 
functions w and w 1 are solutions of the following Cauchy problems 

won = A^ + Att5 (x - x , t) , (x, t) G R 3 x (0, T) , (2.7) 
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w (x,x o ,0) = w ot (x,x o ,0) = 0, (2.8) 

w m = A xWl + a (x) wott (x, x , t) , (x, t) e R 3 x (0, T) , (2.9) 

w (x,x ,0) = w t (x,x o ,0) = 0. (2.10) 
Consider the function h (x, xq, t) , 



t T 

h(x,x ,t) — J dr J Wi(x,x ,s) ds. 





Integrating (2.9) with respect to t twice and using (2.8) and (2.10), we obtain 

h tt = A x h + a(x)wo(x,x ,t),(x,t) 6l 3 x (0,T), (2.11) 

/i(x,a;o,0) = /i t (a;,a;o,0) = 0, (2.12) 

It follows from (2.7), (2.11), (2.12) and the formula (7.13) of §1 of Chapter 7 of [18] that the 
function h(x,xo,t) is 

h(x,x ,t) = — — — I ^- / \y - x \ 2 a(y)dujy, (2.13) 

2n (t z - \x-x \ ) J 

S(x,xo,t) 

where du y = sin9d(pd8, {tp,Q) £ [0, 2n] x [0,7r] are angles in the spherical coordinate system 
with the center at {x } and S(x,x ,t) is the following ellipsoid with foci at {x} and {x } 

S (x, x , t) = {y e R 3 : \x - y\ + \x -y\=t) . 

Setting in (2.13) xo '■— x and denoting v(x,t) = h(x,x,t) , we obtain that the function 
v is the spherical Radon transform of the function a, 

v{x,t) = ^ J a(y)dw y . (2.14) 

\x-y\=t/2 

On the other hand, (2.14) implies that the function v(x,t) = v (x, 2t) ■ t is the solution of 
the following Cauchy problem 

v tt = Av, (x, i)el 3 x (0, 2T) . (2.15) 

v \ t =o= 0,v t \ t =o= a(x) . (2.16) 

Also, using the above linearization one can "translate" the data p(xo,t) in (2.4) for the 
Inverse Problem 3 in the following function p(x,t) 

v\s T =p(x,t),te(0,2T). (2.17) 
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Since the function a (x) — outside of the domain Q, then solving the initial boundary value 
problem (2.15)-(2.17) for (x,t) G (IR 3 \fi) x (0,T) , we obtain the normal derivative q (x,t) , 

^ \s T =q(x,t),te(0,2T) (2.18) 

In conclusion, we have reduced the linearized Inverse Problem 3 to the tp— problem, 
which consists in the recovery of the function a (x) from conditions (2.15)-(2.18). A similar 
derivation is valid for a similar inverse problem for the telegraph equation (2.6) at a = 0, 
see [18] and [21]. 



3 The Method 

We consider Inverse Problem 1, because it is more general than Inverse Problem 2. De- 
note Mu = u tt — Lu. To solve the Inverse Problem 1 numerically, consider the Tikhonov 
regularizing functional 

Je (u) = \\Mu - F\\1 2(Qt) + e |M|^ 2(Qt) 

+ \\Df*u \ St -D^f\\l 2{ST) + ||«, \s T -9\\1 2{St) (3.1) 

+X V IMM) - ^llL(n) +X4, \Hx,0) - ip\\ 2 Hl{n) yu G H 2 (Q t ). 
Here e > is the regularization parameter, 

1 ® u 1 

U u \s T -= 7^ I S T 

and D 13 , \/3\ < 1 is the operator of (x, t) derivatives with, where ^-derivatives are those, which 
are taken in directions orthogonal to the normal vector. Also, 

J 1 for the ip — problem 1 J 1 for the ip — problem 

^ \ for the ip — problem J' X ^ \ for the ip — problem 

Hence, XipX^p — 0, X^ + — 1- I n previous works on the QRM terms in the second line of 
(3.1) were absent because of subtracting off boundary conditions from the original function 
u. Terms in the third line of (3.1) were absent also, and they are incorporated now to 
emphasize the knowledge of one of initial conditions. 

To find the minimizer of J £ (u) , we set the Frechet derivative of this functional to zero 
and obtain for all v G H 2 (Qt) 

J MuMvdxdt + J (D^vD^u + vu) \ St dS + J (v u u u ) \ St dS 
+ Xii> J [V-uVw + uv] (x, 0) dx + Xlp J u t(x, 0)v t (x, 0)dx + e [u,v] (3.2) 
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= J FMvdxdt + J ( D ^ v \s T )D^fdS + + J (v u \s T )-gdS 

Qt S t l^l- 1 S t 

j [Vy?Vv (x,0) + (pv (x,0)} dx + x^J ifr)t(x, fydx. 

n n 
Riesz theorem and (3.2) imply 

Lemma 3.1. For any vector function (F, f, g) G L 2 (Qt) x H 1 (St) x L 2 (St) there exists 
unique solution u e e H 2 (Qt) of the problem (3.2) and 



ll u e|lH 2 (g T ) - ()\ F \\l 2 (q t ) + II/IIhi(s t ) + IMIl 2 (s t ) + X$ II V 9 II m (n) + X v IMIl^q)) ■ 

Setting in (3.2) v := u, we obtain that the unique minimizer of the functional J £ (u) 
satisfies the following estimate 

W Mu Wl 2 (q t ) +Xv \\ u ( x >°)\\hHci)+X<p IK(^,0)llL(n) 

+ \\ U IstIIh1(S t ) + IK \st\\l 2 (S t ) ( 3 - 3 ) 

^ II F IIL(Qt) + II/IIhi(s t ) + \\9\\l 2 (s T ) + ^ IMIffi(n) + X v U\?m m ■ 
To prove convergence of our method, we need to derive from (3.3) the Lipschitz stability 
estimate for the function u in the H 1 (QT)-norm. This in turn requires a modification of 
the proofs of [4], [10], [11] and [13]. We specifically refer to the proofs of Theorem 2.4.1 in 
[13] and Theorem 4.4 in [4]. The main difference with previous proofs is that now either 
\\u(x, 0)\\ 2 H i( n ) or ll' u t( x ; 0)llz, 2 (fi) can b e estimated via the right hand side of (3.3), which was 
not done before. This is because terms in the third line of (3.1) were not included in the 
Tikhonov functional for QRM. So, if x%j> = 0> then x v — 1 an d we estimate \\u(x, 0)\\ 2 L2 ^ in 
the ^—problem. If, however, Xtp — 1> then x^ — and we estimate ||w t (x, 0)|| L2 ^ in the ip 
problem. It is because of the incorporation of these terms that we can assume that T > R, 
as it is required in Inverse Problem 2 (see (1.9)), instead of T > 2R of previous works. The 
required modification is done in the next section. 

4 Lipschitz Stability Estimate 

Theorem 4.1. Let Q C K" be a convex bounded domain with the piecewise smooth boundary 
and let T > R. Suppose that the function u e H 2 (Qt) satisfies the inequality 

W Mu Wl 2 (q t )+X^ \Hx,0)\\ H1(n) + x v ht(x,0)\\ L2{n) (4.1) 

+ \\ u I St 1 1 //i (St) + IK \s t \\l 2 (S t ) - K -> 

where K = const. > 0. Then 

\\ u \\m(Q T )+X v \\u(x,0)\\ H1(n) + x^ IK(a:,0)|| La(n) < CK. (4.2) 



10 



Proof. Choose a pair of points x',y' G <9f2 such that \x' — y'\ = 2R. Put the origin at the 
point (x + y) /2. Choose a constant rj G (0, 1). Consider the function p (x, t) , 

p (x, t) = |x| 2 — r/t 2 . 

Consider the Carleman Weight Function (CWF) C(x, t), 

C(x,t) = exp(2Ap(x,£)) , 

where A > 1 is a parameter. For any positive number b denote 

Gb = {( x ,t) | p (x, t) > b, x G f2, t > 0} . Choose a sufficiently small number c G (0, R 2 ) . 
Since T > R, then in G c 

t 2 <^^<T 2 ,V^G( % ,l), 

where rj = r/ (T, i?) G (0, 1) . Hence, G c C Qt- Choose 5 G (0, c) so small that G c+ as ^ 0- 
Note that 

G c +4<5 C G c +3<5 C G C+ 2S C G c +<5 C G c . (4-3) 

Denote M u = u tt — Au. The following pointwise Carleman estimate takes place 
{M u) 2 C 2 > CX {\V x , t u\ 2 + AV) C + V • U + V t in G c , Vm g C 2 (G c ) , VA > A , (4.4) 



where 



+ \V\ < CX (\V x>t u\ 2 + AV) C in G c 
and Ao (G c , rj) > 1 is sufficiently large. In addition, the function V is estimated as 

\V\ < CXh (| V^mI 2 + u 2 ) C + CA 3 |«t| (I Vu| + |u|) C in G c . 



(4.5) 



(4.6) 



This Carleman estimate was proven in Theorem 2.2.4 of [13]. It was derived earlier in §4 of 
Chapter 4 of [18]. 

Consider the cut-off function p (x, i) G C 2 (G c ) such that 



^ 1 in G c+ 28 

p(x,t) = <{ in G c \G c+5 

between and 1 otherwise 



> . 



(4.7) 



The existence of such functions is well known. For an arbitrary function u G C 2 (G c ) denote 
v = v(u) := pu. Using (4.4)-(4.6), we obtain 

J (M v) 2 Cdxdt>CX J (\V x>t v\ 2 + X 2 v 2 )Cdxdt 



-CX 3 / [\u t \{\Vu\ + \u\)}{x,0)exp(2X\x\ 2 )dx-CX / \{D p u) 2 + X 2 u 2 u 



CdS. 



n 



St 



n 



Because by (4.3) G c+2 8 C G c , then (4.7) implies that the last inequality can be rewritten as 
J (M v) 2 Cdxdt >CX J (\V x , t u\ 2 + \ 2 u 2 )Cdxdt (4.8) 



c+2S 



CX 3 / [\u t \(\Vu\ + \u\)}(x,0)exp(2X\x\ 2 )dx-CX / \(D p uf + u 2 v 



CdS. 



By (4.7) the left hand side of (4.8) can be estimated from the above as 

J (M v) 2 Cdxdt < J (M u) 2 Cdxdt + C J (IV^m) 2 + u 2 ) Cdxdt 



G C+ 2S 



G c \G c+2 g 



< J (M u) 2 Cdxdt + CX 3 |H|^ 1(QT) exp[2A(c + 25)] 

G c +26 

< J {Muf Cdxdt + C J (\V X;t u\ 2 + u 2 ) Cdxdt + C\\u\\ h1(Qt) exp [2X(c + 25)]. 

G c +28 G c+ 25 

Substituting this in (4.8), recalling that A is sufficiently large and using (4.3), we obtain 

J {Muf Cdxdt + A 3 \\u\\ 2 h1(Qt) exp [2A (c + 25)] 



+A 



G C+ 2S 

,2 



(D^u) + A 2 m 2 / J CdS + A 3 J []u t \ (| Vu\ + \u\)\ (x, 0) exp (2A \x\ 2 ) dx 

S T n 

>CX J (\V x , t u\ 2 + X 2 u 2 ) Cdxdt > CA 3 exp [2A (c + 35)] J (\V x , t u\ 2 + u 2 ) dxdt. 

G C +2S G c+ ss 

Let m = maxg c p(x, t). Dividing the last inequality by CA 3 exp [2A (c + 35)] , we obtain 

J (\V x , t u\ 2 + u 2 )dxdt < 



Gc+36 



Ce 2Xm (||M M || 2 2(Qt) + || u \ St \\ 2 HHSt) + IK \sX 2 ( St) ) +C\\u 



1 2 2X8 
\H^{Q T ) e 



(4.9) 



+Ce 2Xm J [\u t ] (| Vu| + \u\)\ (x, 0) dx. 
n 
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The last term of (4.9) was not present in previous publications, and we will analyze it 
now. Consider the ^—problem first. That is, consider the case x v — 1) Xip — 0- Let 7 > 
be a small number which we will choose later. We estimate the last term of (4.9) as 

Ce 2Xm J [\u t \ (\Vu\ + \u\)\ (x, 0) dx 
n 

(\Vu\ 2 + u 2 ) (x, 0) dx + — ^— J u 2 (x, 0) dx (4.10<f) 



n n 

n g4Am 

<C 7 ||«(x,0)||i 1(n) + — —K 



2 



We have used (4.1) to estimate the last term in the second line of (4.10y?). Consider now 
the ip— problem. Then similarly with (4.10</?) 

[\ut\ (|Vu| + \u\)} (x, 0) dx < C 7 |K (x, 0) \\ 2 L2{n) + K 2 . (4.1(ty) 

Q 

Consider now the set 

F 1 (c,5) = G c+35 n{te(0,5)}. 

Then 



{x,t) : \x\ > yc + 35 + r)5 2 ,x e n,t e (0,5) } c F 1 (c,5) . ( 1.11 ) 

Then (4.9), (4.10y?) and (4.10^) imply that 

J {\^xM 2 + u2 ) dxdt < — K 2 + C Mlun^ e~ 2XS 



Ce 4\m 



Fi(c,S) 



■\HHQ T )' 



+Ce 2Xm (\\u \s t \\ 2 h i (St ) + IK \st\\1 2 (s t )) 
+X v C-i \\u (x, 0)\\ 2 H1(Q) + x^Cl IK ( x > °)\\l 2 (Q) ■ 

Hence, by (4.1) 

//nr 4Am 
(IV^m) 2 + u 2 ) dxdt < K 2 + C \\u\\ 2 h1{Qt) e- 2XS (4.12) 
T 

+X V ^7 Ik O)IIhi(O) + ^7 IK fa, °)llL(n) ■ 

Choose numbers c and 5 so small that 3a/c + 35 + < R. Hence, we can choose x £ 
such that |x | = 3a/ c + 35 + ?7<5 2 . Next, we "shift" the function p(x,t) to the point x , thus 
considering the function 

p (x — Xq, t) — \x — xq\ 2 — i]t 2 . 



Fi(c,S) 
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For b > let Gb (x ) = {(x,t) \ p (x — x , t) > b, x E f2, t > 0} . Similarly with the above 
denote 

F 2 (c,5) = G c+35 (x )n{t e (0,5)}. 

Then 

j(x,t) : \x-x \ > ^Jc + 35 + r]5 2 ,x e Q,t E (0,5) | C F 2 (c,S). (4.13) 
Consider an arbitrary point x such that x G < 2a/ c + 35 + r/5 2 ^ . Then 

|x - x| > |x | - \x\ > 3\Jc + 35 + r]5 2 - 2\Jc + 35 + i]5 2 = \Jc + 35 + r)5 2 . 
Hence, by (4.13) 

j(:M) : |x| < 2^c + 35 + r]5 2 ,t £ (0,5) J C F 2 (c, 5) . 

Combining this with (4.11), we see that 

Q x (0, 5) = Q 5 C F 1 (c, 5) U F 2 (c, 5) . (4.14) 
Using function p (x — x , t) instead of p (x, t) , we obtain similarly with (4.12) 

{\V X ,M 2 + u 2 ) dxdt < K 2 + C \\u\\ 2 h1{Qt) e- 2 " 5 



F 2 (c,8) 

+X^7 \\u (x, 0)||^i (n) + x^l \\u t (x, 0)\\ 2 L2{Q) 
Combining this with (4.12) and (4.14), we obtain 

+X v C-f \\u (x, o)\\ 2 m{n) + x*Cl IK (x, o)\\l 2{n) 

Hence, there exists a number t* £ (0, 5) such that 

Ce 4Am „ C 



J (|V X)t u| 2 + « 2 ) (x,t*)dx < 



-I- — \U,\\ £ p-2A5 

S 1 K + S " U " Hl (Qr) e 



1 



X<pCl \\u (x, 0) || Hl(n) + x^C-i \\ut (x, 0) || La(n) 
This inequality combined with (4.1) and the standard energy estimates implies that 

||w|Ihi ( q t ) + X v \\u (x, 0)\\ 2 H1(n) + x^ |K ( x > °)llL(n) ^ C H M I' 2 ' " A<> 
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\hhq t ) e 



Ce 4Xm 

+ —^ R2 + X v C-y W u °)Hffi(n) + Xi,Cl \\u t (x, 0)\\l 2{n) . (4.15) 
Note that 5 is independent on A. Choose sufficiently large A such that 

1 _ Ce - 2X « S > - 
2 

and set A := A . Choose 7 so small that C7 < 1/2. Then we obtain (4.2) from (4.15). □ 

5 Convergence 

Theorem 4.1 enables us to prove convergence of our method. Following the Tikhonov concept 
for ill-posed problems [18], we first introduce an "ideal" exact solution of either <p or ip 
problem without an error in the data. Next, we assume the existence of the error in the 
boundary data / and g and prove that our solution tends to the exact one as the level 
of error in the data tends to zero. We consider the more general Inverse Problem 1. Let 
/* G H 1 (St) and g* G L2 (St) be the exact boundary data (1.3), F* G L2 (Qt) be the exact 
right hand side of equation (1.1) and (p* and ip* be exact initial conditions. We assume that 
there exists an exact function u* G H 2 (Qt) satisfying 

u* tt — L(x, t)u* + F* (x, t) in Q T , (5.1) 

with initial conditions 

u* (x, 0) = if* (x) , u* t (x, 0) = vp* (x) , (p* G H 1 (Q) , V* G L 2 (SI) , (5.2) 

Ou* 

u* \s T = f* ( x , » \s T = 9* (x, t) , (5.3) 

where pf and ip* are exact initial conditions. We assume that the real boundary data in 
(1.3) have an error, so as the given initial condition. In other words, we assume that 

11/ " f*W H Hs T ) + \\9 ~ 9*\\ L2{St) + 11^ - ^*IIl 2Wt) (5-4) 

where 5 > is a small number. The following convergence theorem holds 

Theorem 5.1. Suppose that T > R. Let u £ s G H 2 (Qt) be the solution of the QRM 
problem (3.2), which is guaranteed by Lemma 2.1. Let conditions (5.1)-(5.4) be satisfied. 
Then the following estimate is valid 

\\ u ~ u *Wm{Q T ) + X v \\V ~ <P*\\hhq) + U ~ ^W^n) < C (S + Ve) . 

Proof. Since the functional Jo (u) with the exact data (5.2), (5.3) achieves its minimal 
zero value at u := u*, then the function u* satisfies equation (3.2) with e = and with the 
exact data (5.2), (5.3). Subtracting that equation for u* from equation (3.2) for u := u £ $, 
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denoting w = u £ s — u*, setting in resulting equation v :— w and using (5.4), we obtain 
similarly with (3.3) 

J {Mwf dxdt + \\w( x Mhhq)+X v \\ w t( x , )\\l 2 (n) 
Qt 

+ 11™ \s T f H i {ST ) + \\wu \s t \\1 2(St) <45 2 + e. 
The rest of the proof follows immediately from Theorem 4.1. □ 

6 Numerical Implementation 

In our numerical study we have considered the Inverse Problem 2. To generate the data for 
the inverse problem, we have solved the Cauchy problem 

u tt = Am, (x,t) e M 2 x (0,T) , (6.1) 

u(x, 0) = ip (x) , u t (x, 0) = ip (x) . (6.2) 

In our numerical experiments vp (x) = for the ^—problem, and ip (x) = for the -^—problem. 
Because of (1.5) and the finite speed of propagation, we use in our solution of the forward 
problem zero Dirichlet boundary condition at the boundary of the rectangle (— T, a + T) x 
(— T, a + T) (Figure 1). Hence, we solve initial boundary value problem inside of this rect- 
angle for equation (6.1) with initial conditions (6.2) at zero Dirichlet boundary condition. 
In all our calculations we took a = 1. In tests 1, 2 and 5, which are concerned with the 
Inverse Problem 2, we took T = 3. Hence, condition (1.9) is satisfied. Tests 3 and 4 are 
concerned with the Inverse Problem 1 and we have taken different values of T in these tests. 
The square SQ(a) is SQ(a) = SQ(1) = (0, 1) x (0, 1) , the domain Q in tests 1,2 and 5 is 

Q := (0,4) x (0,4) (6.3) 

and in all tests 

<p(x) =t/j(x) =0 forrr^ SQ(1). (6.4) 

We have solved the Cauchy problem (6.1), (6.2) via finite differences using the uniform 
grid. We set 

u(t k , x ln , x 2m ) ~ u kmn , k = 0, N u n = 0, N x , m = 0, N y , 

tk = kht, X\ n = nh xi , X2m = m h X2 > 

step sizes h xi = h X2 = 0.1, h t = 1/15 and N x = N y = 10, N t = 45. This solution has 
generated the boundary data (1.8). Next, we have introduced noise in these data as 

fn {x\ tj) = f {x\ tj) (1 + 7 iV (tj)) , g n {x\ tj) = g {x\ tj) (1 + (t,)) , (6.5) 
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where (x\tj) is the grid point at the boundary. Here N G (—1,1) is a pseudo random 
variable, which is given by function Math.randomQ in Java and 7 G [0.05, 0.5] is the noise 
level. We have chosen the grid points the same as ones in the finite difference scheme we 
have solved the problem (6.1), (6.2). The presence of the random noise in the date prevents 
us from committing "inverse crime". In (6.5) points x % G Tyr U T 2 t- As to T 3 t U T^t, we 
simply set / = g = on this part of the boundary, because of (1.7). 

To find the minimizer of the functional J £ , we have also used finite differences. We have 
used in (3.1) the finite difference approximations for Mu = u tt — Au and u v \s T and have 
minimized the resulting functional J £ with respect to the vector {uk mn } , which approximates 
values of the function u at grid points. Here J e means the functional J £ , which is expressed 
via the finite differences. The norms ||u xi (x, 0)|| i2 ^ , ||u X2 (x, 0)|| i2 ^ in \\u(x, 0)|| H i^ in 

the vp— problem were calculated via finite differences. As to the term \\D^u \s t — ^ f\\L 2 {s T ) 
in (3.1), we have used only (3 — 0, thus ending up with \\u \s T — f\\\ 2 (s T ) ( m the discrete 
sense). Note that since (3 = 0, our numerical results seem to be stronger than Theorem 4.1 
predicts. The integrals were calculated as 

r h h h Nt-lN y -l N x -l 

/ (U tt - Aufdv « £ £ £ M Ln, 

^ 1 fe=l m=l n=l 

where 

M kmn i^k+l,mn ^^kmn ^k— l,mn) Ay i^k,m+l,n ^^kmn ^k,m—\,n) 

A x {Ukm,n+1 ^^kmn ~\~ ^km,n—l) 
.run -l.ran ,m+l,n ,m—l,n) ^x(Ukm,n+l ~\~ ^fcm,n— l) ^tUkmni 



where 



Also, 



,2 



tit ti\ 

n xi n x 2 



/ / (u(t,X 2 ,xl) - f(t,X 2 )) 2 dx 2 dttt h t h X2 J2J2 H km' 

J J 1 — n ™ — n 







k=0 m=0 



where 



H km tlkmn, h km , 

where n* in the layer number (value of x\) at which the grid function f km is given. 

To minimize the functional J e , we have used the conjugate gradient method. Derivatives 
with respect to variables u kmn where calculated in closed forms, using the following formula 

du kmn 

= kk mm n n, 

° a kmn 
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Number of iterations Number of iterations 



Figure 2: Typical dependence of the functional J = J £ (left) and g = ||VJ e || 2 (right) on 
number of iterations. 

where 5 k ^ is the Kronecker symbol. This formula can be conveniently used to obtain closed 
form expressions for derivatives 

dJ £ (ll) 

Let a be the vector of unknowns of the functional J £ . We start our iterative process 
from a := a = 0. It is well known in the field of ill-posed problems that the number of 
iterations can often be taken as a regularization parameter, and it depends, of course on the 
range of parameters of a problem one considers. We have found that the optimal number 
of iterations for our range of parameters is 300. Thus, in all numerical examples below 300 
iterations of the conjugate gradient method were used, thus ending up with a 300 . Figure [2] 
displays typical dependencies of the functional J e (a^) and the norm of its gradient on the 
iteration number k. 



7 Numerical Results 

In this section we present results of some numerical experiments. We have always used 
e = 10~ 6 . Larger values of e such as 10~ 5 brought lower quality results. In our numerical 
experiments we have imaged both smooth slowly varying functions and the finite difference 
analogue of the 5— function. Let (xik,%2r) G be a fixed grid point. To obtain the finite 
difference analogue of 6 (x\ — xik, xi — X2 r ), we consider the following grid points (xi n , xim) 
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Figure 3: Exact (red) and calculated (black) functions ip in (7.1) without balancing coeffi- 
cients with 5% noise in boundary data. 

and we model the function 5 (x\ n — x\ k) x 2m — X2 r ) as 

3 

Qll X2 ll Xl 

where the multiplier at S nk 5 mr is chosen such that the volume of the pyramid based on 
(xi fc _i,x 2r -i) , (x lk -i,x 2r+ i) , (x lk+ i,x 2r +i) , {xik+ux<2r-x) equals to 1. Hence, the support 
of the function S (x ln — xi k , x 2m — x 2r ) is limited only to the point (xi„, x 2m ). 

We have observed that having equal coefficient at all terms of the functional J e in (3.1) 
does not lead to good reconstruction results. This is because not all the terms of (3.1) 
provide an equal impact in this functional. For example, for the (p— problem with no noise 
in the data for the function 

, x _ J sin(27rxi) sin(27TX2), x G SQ (1) 
| otherwise 

we first got the result displayed in Figure [31 One can observe that the error at the boundary 
is significant. And indeed, the values of two terms in (3.1) after 300 iterations were for this 
case 

J {u tt - Aufdxdt w 10~ 3 , ||u - f\\ L2 (s T ) ~ 10" 2 . 
n T 

Hence, the impact of the boundary term in (3.1) is 10 greater than the impact of the 
|| Mw || ^ \ . To minimize the error at the boundary, we took the balancing coefficient 1000 
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Figure 4: Exact (red) and calculated (black) functions ip in (7.1) with balancing coefficients 
with 5% noise in boundary data. 

at 1000 • \\u — /||x 2 (St) instead of 1 • \\u — /||l 2 (5 t )- The other balancing coefficients equal to 
1. The quality of the resulting image was improved, see Figure SJ Thus, in all our tests with 
the ip— problem we have taken the same balancing coefficients. In the case of the ^—problem 
we have taken 100 ■ x^,\\u(x,0) — <p\\ 2 H i^ and the other balancing coefficients equal to 1. 

Note that Theorems 4.1 and 5.1 remain the same, including their proofs, if balancing 
coefficients are introduced. 

Test 1. The ip— problem. Here ip (x) =0 and the function ip (x) to be reconstructed 
is one in (7.1). In Figures [5] and M represent resulting images with 25% and 50% noise 
respectively. Next, we test our method for the case when the term with % is absent in 
the functional J £ (u) in (3.1). Regardless on the small amount of noise in the data, both 
maximal (1) and minimal (—1) values of the imaged function were missed by about 22% 
in this case, whereas they were not missed in the previous cases with 25% and 50% noise 
when the term with x v was not absent in (3.1). To see this, we display on Figure [7] the 
1-dimensional cross-sections by the straight line {x\ = 0.5} of the correct function (7.1), the 
imaged function with 50% noise of Figure [6] and the imaged function with the absent term 
with Xip and 5% noise. One can observe that the maximal value of the calculated function is 
0.7, while the maximal absolute value of the correct function is 0.9, so as the one of Figure 
[HI Here we have 0.9 instead of 1 only because the points with the absolute value of 1 are not 
the grid points. This emphasizes the importance of the incorporation of the term with Xu>- 
We have observed the same for the ip— problem (images not shown). 

Test 2. The ip— problem. In this case ip (x) = and the function ip (x) to be reconstructed 
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Figure 5: Test 1. Exact (red) and calculated (black) functions <p in (7.1) with 25% noise in 
the boundary data. 




Figure 6: Test 1. Exact (red) and calculated (black) functions <p in (7.1) with 50% noise in 
the boundary data. 
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exact 
50% noise 

5% noise 
no integral 




Figure 7: Test 1. Cross sections of exact (red) and calculated (black, blue) functions if with 
5%, 50% noise, "no integral" means X v =o- O ne can see that the maximal value of the case 
X v =o is 0.7/(— 0.7). The maximal value of the exact function is 0.9 < 1 only because of the 
grid step size. 
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Figure 8: Test 2. Exact (red) and calculated (black) functions ip in (7-2) with 5% noise in 
the boundary data. 

is 

. , v / sin(^ (ari - 0.5) sin(f (x 2 - 0.5)), x e SQ (1) 
[ otherwise. 

Figures [HI 191 and [TUHisplay resulting images of the function (7.2) with 5%, 25% and 50% of 
the noise level in the data respectively. 

Test 3. The if— problem in SQ(1) for T G (0.5diamSQ (1) ,diamSQ (1)) , where 
diamSQ (1) = \/2 is the diameter of the square SQ (1). We have decided to see what kind 
of results can be obtained if the boundary Cauchy data are given on the entire boundary of 
the square SQ (1) in the case when T 6 (0.5diam SQ (1) , diam SQ (1)) . We are especially 
interested in the question about the influence of terms with x v an d Xib- We have used 

T = 0.75 < diam SQ (1) = V2, N x = N y = 20, N t = 3. 

and have reconstructed the function (7.1). Figure [TT] displays the resulting image with 25% 
noise in the case when the term \ v is present in (3.1). This quality of the reconstruction is 
good for such a high noise level. Figure [12] displays the 1-dimensional cross-section of the 
image by the straight line {x\ = 0.5}, as well as the 1-dimensional cross-section of the image 
for the case when the term with \ v is n °t present in (3.1) and 25% noise in the data is in. 
One can observe that the minimal value of (—0.9) is not achieved in the case when the term 
with Xip is n °t present. The calculated minimal value is (—0.7) in this case. 

Test 4. The ip— problem in SQ(1) for T > diamSQ (1). We now test our method for 
the case when the boundary Cauchy data are given at the entire boundary of the rectangle 
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Figure 9: Test 2. Exact (red) and calculated (black) functions ip in (7.2) with 25% noise in 
the boundary data. 




Figure 10: Test 2. Exact (red) and calculated (black) functions ip with 50% noise in the 
boundary data. 
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Figure 11: Test 3. Exact (red) and calculated (black) solutions of the problem tp— in SQ(1) 
with 25% noise in the boundary data for T = 0.75. 
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Figure 12: Test 3. Cross sections of exact (red) and calculated (black, blue) functions <p with 
25% noise in the boundary data for T = 0.75, "no integral" means Xu> = 0- The maximal 
value of the exact function is 0.9 < 1 only because of the grid step size. 
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SQ (1) and T > diam SQ (1). We take 



T = 2, N x = N y = 20, N t = 60. 

The function (7.1) was reconstructed. Figure fT3l displays the resulting image with 25% noise 
and Figure [TH displays the 1- dimensional cross-section of the image by the straight line 
{x\ = 0.5}, as well as the 1-dimensional cross-section of the image for the case when the 
term with x v is not present in (3.1) (with 25% noise). One can observe that both images are 
very close to the correct one. This points towards the fact, which follows from the theory of 
above cited publications and also from Theorem 4.1: the presence of terms with x<p an d Xtp 
is important only when T G (R, 2R) and it is unimportant for T > 2R. 

Test 5. The (p— problem with two 5— functions. We now again consider the Inverse 
Problem 2 with the domain Q as in (6.3) and with ip(x) = 0. The data for the forward 
problem were simulated for the case 

(p (xi,x 2 ) = 5{x 1 - 0.4, x 2 - 0.4) +5{x 1 - 0.7, x 2 - 0.7) (7.3) 

with the above described finite difference analogue of the 6— function. Figure [15] displays 
the resulting image of the function (7.3) for the case of 50% of the noise in the boundary 
data, scatter plot mode was used, squares show exact height. Figure [TBI on the other hand, 
shows the image when the term with x v is absent in (3.1) and only 5% noise in the data is 
present. One can see that the correct height is not reached on Figure [TB| unlike Figure [151 
This again shows the importance of the introduction of terms in the third line of (3.1). 

Very similar results (not shown) were obtained for the ^—problem with exactly the same 
5— functions as ones in (7.3). 

8 Conclusions 

We have considered the inverse problems of the determination of one of initial condition 
in a hyperbolic equation using the lateral Cauchy data. We have presented applications of 
these problems to the thermoacoustic tomography, as well as to linearized inverse acoustic 
and inverse electromagnetic problems. The problems we consider are very close ones with 
the Cauchy problems for hyperbolic equations with the lateral data, and we have actually 
solved the latter numerically in Tests 3 and 4. We have focused on the inverse problem in an 
infinite domain (octant), whereas only finite domains were considered in previous numerical 
studies. Nevertheless, we are able to reduce our inverse problem to one in a finite domain 
due to the finite speed of propagation of waves. Since one initial condition is known, we 
were able to decrease the observation time T by twofold. We have shown numerically that 
it is important to know one of initial conditions if T < diameter (Q) , as it is required by 
stability and uniqueness results. However, if T > diameter (Q) , then both the theory and 
our numerical result of Test 4 show that one does not need to know the initial condition. 

We have proposed a new version of the Quasi-Reversibility method. The main new 
element of this version is the inclusion of the terms characterizing a priori knowledge of 
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Figure 13: Test 4. Exact (red) and calculated (black) solutions of the <p— problem in SQ(1) 
with 25% noise in the boundary data for T = 2. 
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Figure 14: Test 4. Cross sections of exact (red) and calculated (black, blue) functions <p with 
25% noise in the boundary data for T = 0.75, "no integral" means Xu> = 0- The maximal 
value of the exact function is 0.9 < 1 only because of the grid step size. 
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Figure 15: Test 5. Exact (red) and calculated (black) function (p with 50% noise in the 
boundary data. The function x v m (3-1) is present. Scatter plot mode. Squares show 
heights. Correct heights are achieved. 
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Figure 16: Test 5. Exact (red) and calculated (black) function ip with 5% noise in the 
boundary data and = 0. Scatter plot mode. Squares show heights. Correct heights are 
not achieved. 
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one of initial conditions. Two other new elements are incorporation of boundary terms 
in the Tikhonov functional instead of subtracting off boundary conditions and the use of 
finite differences instead of finite elements in the inverse solver. To prove convergence of 
this new version, we have modified the technique of previous works, which is based on 
Carleman estimates. A comprehensive numerical study of the proposed numerical method 
was conducted. This study has demonstrated robustness of our technique with respect up to 
50% random noise in the data, similarly with previous publications [4], [12], [15]. This study 
has also demonstrated that this method is capable to image sharp peaks, which is important 
for the application to thermoacoustic tomography, for example. 
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